Thermodynamics of a trapped unitary Fermi gas 
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Thermodynamic properties of an ultracold Fermi gas in a harmonic trap are calculated within a 
local density approximation, using a conserving many-body formalism for the BCS to BEC crossover 
problem, which has been developed by Haussmann et al. [Phys. Rev. A 75, 023610 (2007)]. We focus 
on the unitary regime near a Feshbach resonance and determine the local density and entropy profiles 
and the global entropy S(E) as a function of the total energy E. Our results are in good agreement 
with both experimental data and previous analytical and numerical results for the thermodynamics 
of the unitary Fermi gas. The value of the Bertsch parameter at T — and the superfluid transition 
temperature, however, differ appreciably. We show that, well in the superfluid regime, removal of 
atoms near the cloud edge enables cooling far below temperatures that have been reached so far. 

PACS numbers: 03.75.Ss, 03.75.Hh, 74.20.Fg 



I. INTRODUCTION 

The BCS to BEC crossover problem of a Fermi gas with 
an adjustable attractive interaction has been investigated 
theoretically for quite some time [H, H U 0, [E IE 0- For 
low temperatures the gas is superfluid and, in the case of 
s-wave interactions, it exhibits a smooth crossover from 
the well known BCS regime of weakly bound Cooper 
pairs to the BEC regime of tightly bound bosonic dimers 
with a residual repulsive interaction d [E 0. In re- 
cent years, this crossover has been realized experimen- 
tally using ultracold Fermi gases in optical traps, where 
the interaction can be tuned using Feshbach resonances 
H EE HO, El El EI El EE E3- In the experimen- 
tally relevant case of so called broad Feshbach resonances, 
which is in principle always realized in the dilute limit, 
the physical properties of the homogeneous gas at equal 
densities for both spin components are described by only 
two dimcnsionlcss parameters: the interaction strength 
V = 1/fcpa and the temperature 9 — ksT/sp- Here, 
a is the s-wave scattering length which fully charac- 
terizes interactions in the dilute, ultracold limit, while 
the scales for length and energy are determined by the 
Fermi wave number kp = (S^n) 1 ^ and the Fermi en- 
ergy £p = h 2 kp/2m, respectively, where n = N/V is the 
particle density. 

A particularly interesting regime is located near the 
Feshbach resonance, where the scattering length a is in- 
finite. At this point and, more generally, in the so-called 
unitary regime where kp\a\ ^> 1, the dimcnsionless in- 
teraction parameter v disappears from the problem. All 
thermodynamic quantities are therefore universal func- 
tions of the dimensionless temperature 6 = ksT/ep [lg| |. 
On a microscopic level, the unitary gas exhibits a partic- 
ular kind of scale invariance, similar to a gas with purely 
inverse square two-particle interactions 11911. More gen- 
erally, as shown by Nikolic and Sachdev [20], universal- 
ity is not restricted to the unitary regime. It is tied to 



the fact that the unitary balanced gas at zero density is 
an unstable fixpoint with only three relevant perturba- 
tions. Since there is no small expansion parameter, the 
unitary regime is the most challenging one from a theo- 
retical point of view. In addition, it is in fact precisely 
this regime which is accessible experimentally (see e.g . 
the recent review articles by Ketterle and Zwierlein [2l| . 
by Bloch et al. [22j and by Giorgini et al. p3j). 

In a recent paper [24 |. we have presented a field the- 
oretic approach for the thermodynamics of the BCS to 
BEC crossover, which is based on the formalism devel- 
oped by Luttinger-Ward [25[ and DcDominicis-Martin 
[261 ]. In the following, this approach for the homogeneous 
gas is applied to calculate the thermodynamic properties 
of the trapped Fermi gas, using a local density approxi- 
mation. We compare our results with a recent experiment 
by Luo et al. [27| and with recent theories [H [22, HE] ■ In 
particular, we provide results for the entropy as a func- 
tion of temperature, which allows to do reliable thermom- 
etry for the trapped unitary gas and also gives a precise 
value for the critical temperature and the associated en- 
tropy per particle. In addition, we show that starting 
well in the superfluid regime, much lower temperatures 
and entropies can be reached by removing atoms from 
the edge of the cloud, which carry most of the entropy. 



II. LOCAL DENSITY APPROXIMATION 

In our previous paper we have calculated the thermo- 
dynamic quantities for the homogeneous system. At a 
given particle density n = N/V, these are the internal 
energy per particle u = U/N and the entropy per parti- 
cle s = S/N. The Fermi wave number kF = (S^n) 1 / 3 
and the Fermi energy ep = h 2 k F /2m can be used as scale 
factors in order to make the thermodynamic quantities 
dimcnsionlcss. 

In the presence of an external confining potential V(r), 
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the particle density n(r) is non- uniform. Within a lo- 
cal density approximation (LDA), thermodynamic quan- 
tities like pressure p(r), chemical potential /i(r), entropy 
per particle s(r) or the internal energy per particle u(r) 
are then also spatially varying, being determined by the 
corresponding equilibrium values in the uniform system 
evaluated at the local density n(r). The local density ap- 
proximation neglects the dependence of thermodynamic 
properties on density gradients. At zero temperature, it 
is essentially a zeroth order semiclassical approximation 
[31I ]. It is valid as long as the local Fermi wave num- 
ber kp(r) times the oscillator length £q = (h/muj) 1 ^ 2 
defined by the characteristic frequency lu of the confin- 
ing potential is much larger than one. Except near the 
edge of the cloud, where the density approaches zero, 
this condition is well justified for most of the experi- 
ments because typical Fermi energies e f are of the order 
of several kHz, while the trapping frequencies are around 
lu f=a 100 Hz or smaller. Specifically, near the trap center 
kp(0)£o — iV 1 / 6 , where N is the total particle number in 
a trap. As will be shown below, the finite size corrections 
to the ground state energy are of relative order (3N)~ 2 / 3 
in a harmonic trap. They arc therefore negligible, at least 
for global obscrvables, for the typical particle numbers in 
experiment, where N w 1.3 x 10 5 [27[. This conclusion is 
also supported by a recent comparison of LDA with a nu- 
merical solution of the Bogoliubov-DeGennes equations 
& 

Within the LDA, the global thermodynamic quantities 
of the trapped Fermi gas are obtained by integrating over 
the whole space. Specifically, we define 
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N = I d 3 r n(r) , 
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d 3 r n(r) V(r) , (2.2) 
d 3 r n(r) u(r) , (2.3) 
U + E pot = J d 3 r n(r) [u(r) + V(r)] , (2.4) 



d 3 r n(r) s(r) , 



(2.5) 



as the particle number, the potential energy, the internal 
energy, the total energy, and the total entropy. Here the 
particle density acts like a distribution function to define 
averages over the trap. 

In an optical trap where the laser intensity profiles are 
Gaussian functions, the confining potential V(r) is given 
by an anisotropic Gaussian function. Close to the center 
of the trap this potential can be approximated by an 
anisotropic harmonic function 

V(r) = l -m(iolx 2 + toy + lu 2 z z 2 ) (2.6) 

where m is the mass of the atoms and u) x , ui y , lu z are the 
harmonic oscillator frequencies in the three spatial direc- 
tions. In the following, we use the harmonic potential 



(12. 6|) and neglect the anharmonic terms. Since in LDA 
/i(r) is the chemical potential relative to the potential 
V(r), in thermal equilibrium the total chemical potential 
A*tot = A*( r )+^( r ) is constant which implies the condition 



H(t) + V(r) = p(0) 



(2.7) 



This equation together with the requirement of a con- 
stant temperature T determines the spatial dependence 
of all local thermodynamic quantities. The particle den- 
sity n(r) implies a local Fermi wave number kp(r) and 
a local Fermi energy ef(t). As a consequence the di- 
mensionless parameters v(r) = \/kp(v)a and 8(r) = 
ksT / £f{t) depend on the local position in the trap. 

It is convenient to define the weighted radial coordinate 
r by 
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(2.8) 



where lu = (ojxLUyLUz) 1 / 3 is the average harmonic fre- 
quency. With this definition the confining potential ac- 
quires the simple form V(r) = V(r) = ^mu} 2 r 2 . As a 
consequence, the anisotropy of the trap becomes irrele- 
vant because all local quantities n(r) = n(r) etc. only 
depend on the weighted radial coordinate r. Rewriting 
the space integrals in Eqs. (|2.1[) - (|2.5p in terms of r, this 
applies also to all thermodynamic quantities derived from 
the local density. A convenient measure for the overall 
length scale in a trap is provided by the Thomas-Fermi 
radius Rtf = (24A^) 1 / 6 (/i/mw) 1 / 2 of the confined non- 
interacting Fermi gas at zero temperature with a given 
total particle number N . Similarly, as a characteristic 
scale for the energy we define the corresponding Fermi 
energy Ep = (3Ny^ 3 hLu of the non- interacting gas. 



III. UNITARY REGIME AT ZERO 
TEMPERATURE 

The thermodynamic quantities can be made dimcn- 
sionless by considering the ratios /u(r)/ep(r), ii(r)/ejr(r), 
and s(r)/ks- These ratios depend on the space coordi- 
nate r only implicitly via the dimensionless parameters 
v(r) = l/fcp(r)a and 9(r) = fe^T/e^r). A particular 
case is the unitary gas at zero temperature, where both 
parameters vanish identically v(r) = 0, 8(r) = 0. Using 
standard thermodynamic relations, it is straightforward 
to show that all thermodynamic quantities can be ex- 
pressed in terms of a single dimensionless parameter, the 
so-called Bertsch parameter £ [HI, HH in the form 



M(r) 
e F (r) 



£ F {r) 



(3.1) 



These relations hold both for an ideal Fermi gas, where 
£ = 1 and also at the unitarity point, where £ has a 
nontrivial value [TH, [2(| ■ In our previous work [24[ , we 
have calculated both ratios independently and obtained 
n(r)/e F (r) = 0.358 and u(r)/e F (r) = 0.210. The first 
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ratio implies a Bertsch parameter £ = 0.358, while the 
second ratio implies £ = 0.351. These two values differ by 
about 2.0%. Consequently, the relation u(r)//z(r) = 3/5, 
which is valid both for an ideal and a unitary Fermi gas, 
is satisfied up to an error of 2.0%. 

Inserting the harmonic potential (|2.6|) into Eq. (|2.7j) 
and using the weighted radial coordinate (|2.8p , we obtain 
the chemical potential 



M (r) = M (0) [1 - r 2 /r, 



TF\ 



(3.2) 



where rpp is the Thomas Fermi radius of the unitary 
gas at zero temperature. The dimensionless ratios (|3.1[) 
imply similar functional forms for the other quantities, 
e.g. fc F (r) = fc F (0) [1 — r 2 /r^ F \ 1 / 2 for the local Fermi 
wavevector. These expressions are valid only for r < ttf, 
because the particle density n(r) is nonzero only in this 
case and zero otherwise. Evidently, the Thomas-Fermi 
radius rpp of the unitary Fermi gas is the only parameter 
which determines the spatial dependence of the thermo- 
dynamic quantities. It is related to the Thomas- Fermi 
radius Rp f of the ideal Fermi gas by the Bertsch param- 
eter via 



Ttf/R 



TF 



[fx(0)/s F (0)] 



1/4 



1/4 



(3.3) 



Using the dimensionless ratio /x(r) /e F (r) = 0.358 of 
our numerical calculations [24[ we obtain the result 
ttf/Rtf = 0.773. This value agrees very well with 
the most recent field theoretic result £ = 0.367(9) for 
the Bertsch parameter, obtained from a Borcl rcsumma- 
tion of an expansion around the upper critical dimension 
four, carried out to three loop order [35|, It is some- 
what smaller, however, than the result ttf/Rtf = 0.80 
found by using the Bertsch parameter £ = 0.42(1) that 
follows from variational Monte Carlo calculations [53, [38| 
or from a theory that includes the Gaussian fluctuations 
around the BCS mean field, extended to arbitrary cou- 
pling, where £ = 0.40 [H,^. 

Next we insert the local internal energy per particle 
u(r), the potential V(r) = \mJ 1 r 1 , and the local particle 
density n(r) into Eqs. (|2.2p - (|2.4j) in order to calculate the 
energies of the unitary Fermi gas in the harmonic trap. 
We thus obtain the dimensionless ratios 
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£f(0) 
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1/2 



m(o) 



£f(0) 



-1/2 



U + E, 



pot 



NEp 



= 0.444 . 



(3.4) 

(3.5) 
(3.6) 



The explicit numbers are obtained by inserting the ra- 
tios fi(r)/e F (T) = 0.358 and u(r)/e F (r) = 0.210 of our 
numerical calculation [24[. 

For a harmonic potential V(r) it is well known that 
the internal energy U and the potential energy E pot are 



related to each other by the virial theorem U — E pot . 
As shown by Thomas et al. [4(j, this theorem also holds 
for the interacting Fermi gas in the unitary regime. As 
a consequence, the results of Eqs. (|3.4p - (|3.6[) should be 
equal. By using the dimensionless ratios of the homoge- 
neous system (|3.1[) we can express the results in terms of 
the Bertsch parameter £ according to 



E 



= 2- 



U 



NE F NE F 



= 2 



E 



pot 



NE F 
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(3.7) 



In practice, our resulting numbers differ by about 2.0%. 
This difference is related to the fact, that the exact rela- 
tion u(r)//j,(r) = 3/5 for the unitary gas is satisfied only 
with an error of 2.0% in our theory. 

The result (|3.7p for the ground state energy in the trap 
is based on using LDA and provides the exact leading or- 
der contribution in the limit N — > oo. The question of 
how large the subleading corrections to this result are 
has been addressed by Son and Wingate [il[. Using a 
gradient expansion of the effective field theory describ- 
ing the low energy physics of the fermionic superfluid at 
unitarity, they have determined the g 2 -corrections to the 
density response, which is equal to the uniform compress- 
ibility dn/dn at q = 0. These corrections give rise to an 
additional contribution 1411 



E 



NE F 



1 + in- 
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;C 2 - Ci 



u)± + ujf. + lu; 



(3JV) 



-2/3 



(3.8) 



to the ground state energy, which contains the two di- 
mensionless coefficients c\ and C2 which appear beyond 
the leading coefficient £ of the uniform system in an ex- 
pansion up to second order in gradients. The leading 
correction to LDA is thus of relative order (37V) ~ 2 / 3 
(formally it is ~ h 2 ) and describes a curvature instead 
of the naively expected surface correction, which would 
scale like TV -1 / 3 . The absence of a surface correction 
also apppears for an ideal Fermi gas and is a peculiarity 
of the harmonic confinement [3l|. The prefactor of the 
(3N)~ 2 / 3 correction in (|3 . 8|) has been determined from an 
expansion around dimension four by Rupak and Schafer 
[42j and is 2.41 for an isotropic trap. For experimentally 
relevant particle numbers N ~ 1.3 x 10 5 , the beyond LDA 
relative corrections to the ground-state energy are there- 
fore only around 4.5 x 10 -4 and thus are clearly below 
the present experimental accuracy. 

The final results of this section are the Thomas-Fermi 
radius rxF /Rtf = 0.773 and the total energy E/NEp = 
0.444 of the unitary Fermi gas at zero temperature in a 
harmonic trap, which depend only on the Bertsch param- 
eter £ = 0.358 of the homogeneous system by Eqs. (|3.3[) 
and (|3.7p . These results will be compared with experi- 
ments and other theories in the next section, where we 
discuss the situation at finite temperatures. 
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IV. UNITARY REGIME FOR NONZERO 
TEMPERATURES 

For nonzero temperatures, the local thermodynamic 
quantities are space dependent via the scale factors 
fci?(r), £p(r), and also via the dimensionlcss temperature 
9(r) = fcsT/eF(r). In the unitary regime the dimension- 
less interaction v(r) = I/Uf (r)a — is constant. Since 
the effective temperature increases towards the edge of 
the cloud, the local particle density n(r) does not follow 
the simple [1 — r 2 /r^ F ] 3 / 2 law valid at T = and Eq. 
(|2.7|) must be solved numerically in order to obtain the 
detailed density profile n(r) = n(r) as a function of the 
weighted radial coordinate r. 

Using the results of our previous numerical calculations 
[HI for the homogeneous system, we obtain the density 
profiles n(r) for several values of the temperature T which 
are shown in Fig. [1] The black solid line is the density 
profile for zero temperature. It is nearly identical to the 
expected zero temperature Thomas-Fermi profile. Differ- 
ences are due to the limited numerical accuracy, which 
are, however, much smaller than the thickness of the line. 
At zero temperature, the Fermi gas is superfluid in the 
whole trap. The colored and dot-dashed lines represent 
the density profile for nonzero temperatures. For very 
low temperatures, changes occur only close to the sur- 
face of the cloud where v/Rtf ~ ttfIRtf = 0.77. For 
this reason, the brown short dashed line is visible only 
for 0.65 < t/Rtf S 0-85 while otherwise it is nearly on 
top of the black solid line. The red long dashed line is 
still very close to the black solid line. At a position in 
space r the Fermi gas may be normal fluid or superfluid 
if the local dimensionlcss temperature 9(r) = fc^T / 'ep(r) 




FIG. 1: (Color) The local particle density n(r) of the unitary 
Fermi gas in a harmonic trap as a function of the radius r 
for several temperatures 6(0) = ksT '/ 'ef(0) = 0.000 (black 
solid), 0.032 (brown short dashed), 0.065 (red long dashed), 
0.100 (orange short dashed), 0.141 (green long dashed), 0.192 
(blue dot-dashed), 0.252 (magenta doubledot-dashed), and 
0.313 (turquoise dot-doubledashed) . The length scale is 
Rtf = (24N) 1/6 (h/mco) 1/2 . 



is above or below the critical value 9 C s» 0.16. For the 
brown, red, orange, and green dashed lines there exists a 
respective weighted radius r c so that 8(r c ) = 9 C . In these 
cases the Fermi gas is superfluid in the inner region of the 
trap where r < r c and normal fluid in the outer region 
where r > r c . For the blue, magenta, and turquoise dot- 
dashed lines the dimensionlcss temperature is 9(r) > 9 C 
for all positions in space r so that the Fermi gas is normal 
fluid in the whole trap. 

Since 9(r) has its lowest value for r = 0, the dimension- 
less temperature at the center of the trap 9(0) determines 
the superfluid transition of the confined Fermi gas. For 
9(0) > 9 C the Fermi gas is completely normal fluid, while 
for 9(0) < 9 C there exists a superfluid region close to the 
center. 

In a homogeneous gas, the normal to superfluid tran- 
sition is a continuous phase transition of the 3D XY type 
along the complete BCS to BEC crossover, because the 
broken symmetry is that associated with a complex scalar 
order parameter. By contrast, our approach [24| pre- 
dicts a weak first-order superfluid transition because the 
superfluid phase of the Luttinger-Ward theory does not 
smoothly connect with the normal-fluid phase at a single 
critical temperature 9 C = ksT c /eF- As a result, there 
arc two slightly different critical temperatures # CiU pper 
and 9 c, lower- The upper value 9 CAlppcr is defined by the 
condition that for temperatures above this value the su- 
perfluid order parameter vanishes, while the lower value 
#c,iower is defined by the temperature below which the 
normal-fluid phase is no longer stable. Fortunately, the 
difference between both temperatures, which should van- 
ish in an exact theory, is rather small over essentially the 
whole BCS to BEC crossover. In particular, at unitarity, 
the upper and lower values for 9 C are 0.1604 and 01506, 
which is within the present numerical uncertainties in the 
determination of the critical temperature of the unitary 
gas. Indeed, our critical temperature agrees very well 
with the most precise calculations of 9 C so far by Quan- 
tum Monte Carlo calculations, which give 9 C — 0.152(7) 
for the uniform gas at unitarity [43[, a value that has 
been confirmed very recently [44J . 

The existence of two different critical temperatures 
leads to a multivaluedncss of thermodynamic quantities, 
which is an artifact of the first order nature of the su- 
perfluid transition within our theory. In order to avoid 
multivalued local density or entropy profiles, we have con- 
nected the normal and superfluid branches with a kink 
at the point, where they are closest, thus providing an 
optimal approximation to the exact continuous profiles 
in a theory which properly accounts for the continuous 
nature of the transition in the infinite system. This point 
is related to the upper value of the dimensionlcss criti- 
cal temperature 8 C , upper- Hence, in Fig. [T] the brown, 
red, orange, and green dashed lines have kinks located 
at r c , uppcr /i? TF = 0.684, 0.577, 0.433, 0.209, respec- 
tively. The inner, superfluid branches of these lines show 
a bulge in the center of the trap, which is appreciable, 
in particular for the green curve. The formation of a 
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FIG. 2: (Color) The local entropy density s(r)n(r) of the 
unitary Fermi gas in a harmonic trap as a function of the 
radius r for several temperatures. The curves are related to 
those in Fig. [T] An additional curve is included for 8(0) — 
k B T/e F (0) = 0.0065 (black short dashed). 



small bulge upon condensation, which first appears near 
the trap center, has indeed been observed experimentally 
[2l| . although the effect in axially integrated profiles is 
rather small. 

Previous theoretical results for the density profile have 
been obtained by Bulgac et al. [3(| within a Monte-Carlo 
approach. For zero temperature, their result is very close 
to our black solid line in Fig.[T] They obtain the Thomas- 
Fermi radius ttf/Rtf = 0.81 (45[. This is larger than 
our value 0.77 because of the larger value £ = 0.43 of the 
Bertsch parameter. For nonzero temperatures (the same 
as those in Fig. []}, Bulgac et al. obtain single valued 
density profiles which agree qualitatively with our results 
but differ somewhat quantitatively. They also observe a 
supcrfluid bulge which, however, is smaller. 

A very interesting quantity is the local entropy per 
particle s[r). From our numerical calculations [24j of the 
homogeneous system, and within LDA, the resulting pro- 
files s(r)n(r) of the local entropy per volume are shown in 
Fig- El for several values of the temperature T. The colors 
of the curves are related to those in Fig.Q] For high tem- 
peratures where the Fermi gas is completely normal fluid 
(blue, magenta, and turquoise dot-dashed curves), the 
entropy density is distributed over the whole trap with a 
maximum at the center. For low temperatures (orange, 
red, and brown dashed curves) the center part is super- 
fluid but the outer part is normal fluid. Consequently, in 
these cases the entropy density is minimum in the center 
of the trap but maximum at a nonzero radius. For very 
low temperatures (black dashed line) , the main contribu- 
tion of the entropy is located close to the surface of the 
atom cloud. The width of this surface layer decreases for 
decreasing temperature and eventually shrinks to zero in 
the zero temperature limit (black solid line). 

In Fig. [2] we have again eliminated multivalued regions 
which are an artifact of the first-order nature of the tran- 
sition in our theory for the homogeneous system by us- 



ing the same criterion as in the density profiles of Fig. [T] 
Indeed, since the normal to superfluid transition is not 
associated with a latent heat, the local entropy density 
will be single valued. An interesting observation that is 
evident from Fig. [2 is that sufficiently below the onset of 
superfluidity in a trapped Fermi gas, most of the entropy 
is located in the normal region near the cloud edge. This 
observation suggests an efficient way to lower the entropy 
further by removing atoms in the boundary layer and 
simultaneously readjusting the trapping potential such 
that the now smaller system has a radius close to that 
beyond which atoms have been removed. This idea is - of 
course - in the same spirit than the standard evaporation 
cooling in the normal state [46[ , however it is much more 
efficient. Indeed, consider starting with a dimcnsionlcss 
temperature 6(0) = 0.065 (red long dashed curve) or - 
cquivalently - an entropy per particle S/Nks = 0.433 
which is close to that reached in current experiments. 
Removing about 42 percent of the particles in the shell 
beyond r = 0.5 Rtf? will lower the entropy per particle 
by a factor 8.8 to S/Nks = 0.049 and the dimension- 
less temperature by a factor 2.6 to 6(0) = 0.025. In 
turn, for an initial temperature 6(0) = 0.192 (blue dot- 
dashed curve) where the cloud is a normal gas, remov- 
ing the same amount of particles will reduce the entropy 
only by a factor 1.6. Removing atoms in the outer shell 
repeatedly thus provides an effective tool to reach very 
low entropies and temperatures. In practice, for atoms 
confined in an optical dipole trap, this may be achieved 
by lowering the depth of the optical trap, as was done 
e.g. previously in experiments realizing a condensate of 
fermionic dimers [Til ] . 

From the local profiles n(r), u(r), and s(r), it is 
straightforward to calculate the the total Energy E and 
the global entropy S by evaluating the integrals in Eqs. 
(|2.4[) and (|2.5[) . In particular, we may eliminate the lo- 
cal value 6(0) of the dimensionless temperature in the 
trap center, to obtain the function S = S(E). Our re- 
sult for the ultracold Fermi gas in the unitary regime 
is shown in Fig. [3] as bluc-grcen-red solid line. To dis- 
tinguish the superfluid and the normal-fluid parts of the 
curve, the solid line is shown in blue or red color, respec- 
tively. Apparently, this curve is continuous and there 
is no particular feature, which indicates the superfluid 
transition. In fact, since dS/dE is just the inverse tem- 
perature, this behavior is expected not only in a trap, 
where none of the thermodynamic functions exhibits a 
singularity, but even in the homogeneous gas, because the 
supcrfluid transition is continuous. As mentioned before, 
however, our theory predicts a weak first-order transi- 
tion. The solid line is therefore multivalued in the inter- 
vals 0.656 < E/NE F < 0.677 and 1.56 < S/Nk B < 1.66, 
which is indicated by the green section of the line. Within 
these intervals the superfluid transition is located. In 
practice, evidently, the multivaluedness is so tiny that it 
can hardly be seen in Fig. [3J In the multivalued region 
the blue, green, and red branch of the solid line therefore 
lie on top of teach other. 
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FIG. 3: (Color) The entropy S as a function of the to- 
tal energy E for the unitary Fermi gas in a harmonic trap: 
present theory (blue-green-red solid), NSR theory of Hu et al. 
[28l . l29l . |47| | (magenta dot-dashed) Monte-Carlo simulation of 
Bulgac et al. [30| (orange dashed), and experimental data of 
Luo et al. [53] (green data points). 



Bulgac et al. have obtained S(E) from their Monte- 
Carlo calculation (30j . which, in Fig. [31 is shown as or- 
ange dashed line. Clearly, the agreement with our the- 
ory is nearly perfect for high temperatures well above 
the supcrfiuid transition. However, for lower tempera- 
tures close to and below the transition, the results differ. 
The deviations are largest for zero temperature where 
5 = 0. They are essentially due to the fact that the 
ground-state energy Eq in both approaches differ. In- 
deed, from Eq. we obtain E a /NE F = 0.444 while 
the corresponding value of Bulgac et al. Eq/NEf = 0.50 
is larger because of the larger value £ = 0.43 of the 
Bertsch parameter. Apart from the deviations in the 
limit of zero temperature, the theories also differ in their 
predictions of the behavior near the superfluid transi- 
tion temperature. In particular, Bulgac et al. [30j | have 
calculated the critical values E c /NEp = 0.50 + 0.32 = 
0.82 and S c /Nk B = 2.15 at the superfluid transition. 
These results are considerably larger than our predictions 
0.656 < E C /NE F < 0.677 and 1.56 < S c /Nk B < 1.66, 
whose uncertainty is due to the multivaluedncss near 
T c . As a result, the value of the critical temperature 
k B T c /Ep = 0.27 of the unitary Fermi gas in a trap in 
the theory of Ref. [3^ is considerably larger than our pre- 
diction 0.207 < k B T c /E F < 0.220. 

Now, as pointed out above, our result for the critical 
temperature of the uniform gas agrees rather well with 
the most precise numerical calculations of this quantity 
[43l , |44| . In addition, it is also consistent with the re- 
cent calculations of Bulgac et al. [48| , which indicate that 
6c = k B T c /ep is less or equal to 0.15(1), again in very 
good agreement with our numbers. This underlines that 
our self-consistent, conserving theory of the BCS-BEC 
crossover (2~H provides a quantitatively reliable descrip- 
tion of the thermodynamics of a balanced Fermi gas near 



unitarity, despite the problems with the first order na- 
ture of the transition in our approach and the absence of 
a small expansion parameter. Within LDA, which pro- 
vides a rather accurate description in the relevant regime 
of particle numbers N sa 10 5 , our critical temperature 
k B T c = 0.21(1) Ep for the trapped gas is therefore ex- 
pected to be close to the exact result. This is consistent 
with a very recent analysis of the experimental data of 
Ref. H3, which indicates a critical temperature very close 
to this value [49[. Previous, much higher values of the 
critical temperature for both the homogeneous or the 
trapped gas that were obtained in different extensions 
of the theory by Nozieres and Schmitt-Rink [2|| [5(j and 
also in numerical calculations [5l], H3| are clearly ruled 
out. 

An alternative approach to the BCS to BEC crossover 
problem has been developed by Hu et al. [28| for the ho- 
mogeneous system and applied to the ultracold Fermi gas 
m a trap [M, S3- This theory is an extension of the ap- 
proach by Nozieres and Schmitt-Rink (NSR) [|[ to the 
superfluid region at low temperatures. While the order 
parameter A is determined by the standard mean-field 
gap equation, the chemical potential is calculated by a 
particle-density equation, which includes condensed and 
noncondensed bound pairs. The extended NSR approach 
is in fact a limiting case of our theory in which the full, 
sclf-consistcntly determined Green functions are replaced 
by their zeroth order form obtained in BCS theory. Hu 
et al. [2^, [4?| have calculated S = S(E) and obtain a 
result which is shown as magenta dot-dashed line in Fig. 
[3j Again, the agreement is very good for high tempera- 
tures in the normal fluid region while, however, for low 
temperatures in the superfluid region there are devia- 
tions. Hu et al. [H, obtain the ground-state energy 
Eq/NEf = 0.47 which is related to the Bertsch parame- 
ter £ = 0.40. Moreover, similar to the theory of Bulgac 
et al. [3(il ] their critical temperature k B T c = 0.25 Ep and 
entropy S c ~ 2.2 Nk B are considerably larger than our 
values. 

More recently, Hu et al. [47| have published a self- 
consistent result for the entropy S(E) which they call 
"GG approximation" . This method is equivalent to our 
approach for the homogeneous system [24[ and gives re- 
sults that are nearly identical to the prediction of our 
present theory including the critical temperature k B T c = 
0.21 E F . 

While in the experimental setup the correct trap po- 
tential is Gaussian, all the curves shown in Fig. [3] are 
calculated for the harmonic potential (|2.6| . Slight differ- 
ences would occur in the high-tcmperaturc regime where 
the energies E and the entropies S are large. However, it 
is important that all curves are calculated for the same 
potential so that they all converge to a single line in the 
high-temperature limit. 

The entropy versus the total energy has been measured 
experimentally for ultracold 6 Li atoms in an optical trap 
by Luo et al. [27] ]. The data are shown as green points 
with horizontal and vertical error bars. Clearly the data 
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agree with all theories in the high-temperature regime. 
The slightly larger entropy values of the experimental 
data for large energies may be due to the fact that the ex- 
periment is performed for a Gaussian potential while the 
theoretical curves are calculated for a harmonic poten- 
tial. Apparently, for low temperatures and low entropies 
the experimental data agree better with the theories of 
Bulgac et al. [13 and of Hu et al. [H, HzJ than with our 
theory. In particular, an extrapolation to zero entropy 
gives a ground-state energy E /NE F = 0.53, consider- 
ably higher than our value 0.444. It is difficult, however, 
to quantify the error involved in such an extrapolation 
because the determination of the entropy from a com- 
parison with its ideal Fermi gas limit reached after an 
adiabatic ramp to magnetic fields B = 1200 G far on the 
BCS side of the crossover becomes increasingly difficult as 
S approaches zero. As pointed out above, a quite sensi- 
tive parameter which distinguishes previous theories from 
our present one is the critical temperature of the super- 
fluid transition and the associated value of the entropy 
and energy Unfortunately, it is difficult to determine 
the critical temperature and thus also the corresponding 
value of the entropy from measurements of S(E). Based 
on the excellent agreement of our value for T c and the 
Bertsch parameter £ with the most precise numerical or 
field-theoretical results and the fact that all thermody- 
namic relations are obeyed at the few percent level, it is 
likely that our present theory gives a reliable description 
of the thermodynamics of the trapped unitary gas, de- 
spite the fact that, superficially, the agreement with the 
data shown in Fig. [3] is not as good as those of previous 
theories. 



V. THERMOMETRY 

The most important parameter for thermodynamic 
properties is - of course - the temperature, which unfor- 
tunately cannot be measured directly. In principle, this 
is possible from the density profile n(r) which changes as 
a function of temperature. For balanced gases, however, 
this method is not very reliable. Indeed, for low tempera- 
tures, the density converges to the simple Thomas-Fermi 
profile. As shown in Fig. [TJ the deviations from such 
a profile at finite temperatures are extremely small for 
the inner part of the atom cloud (see the red and brown 
low-temperature curves, which are nearly identical to the 
black zero-temperature curve). Only close to the surface 
r/Rpp ~ 0.77, small variations with the temperature 
are observed. Therefore, the signal to noise ratio will be 
small in the interesting regime below T/Tp w 0.1. 

A quantity which is much more sensitive to temper- 
ature is the entropy density, shown in Fig. [5] At low 
temperatures, it is peaked near the surface of the atom 
cloud. Unfortunately, the local entropy density s(r)n(r) 
is not accessible experimentally. We therefore consider 
the total entropy S, which has been measured by [27|, as 
discussed above. 




FIG. 4: (Color) The entropy S as a function of the tempera- 
ture T for the unitary Fermi gas in a harmonic trap. The blue, 
red, and green sections represent the superfluid, normal-fluid 
and multivalued branches of the curve, respectively. 



In Fig. [4l we plot the total entropy S in a trap as a 
function for the temperature T in units of the Fermi tem- 
perature Tp = Ep/ftB = {iNY^huj /ks of the trapped 
non-interacting gas. As in Fig. [3J we indicate the super- 
fluid and the normal-fluid region by blue and red color of 
the curve, respectively. Again, a tiny multivalued region 
is observed close to the superfluid transition which is in- 
dicated by the green section of the line. However, the 
three branches in the multivalued region are nearly on 
top of each other. Consequently, the entropy S = S(T) 
is effectively continuous at the superfluid transition as 
expected. Using the curve in Fig. 31 the measurements 
of the entropy by Luo et al. [27| therefore allow to reliably 
infer the related temperatures T. 

Alternatively, we may consider the excess internal en- 
ergy density [u(r) — uo(r)]n(r) where uo(r) is the ground- 
state internal energy per particle at zero temperature. 
The related profiles are similar to those shown in Fig. [5] 
For low temperatures the excess internal energy density is 
peaked close to the surface of the atom cloud. Integrating 
over the space we obtain the excess internal energy U—Uq 
which is related to the excess total energy by the virial 
theorem E — Eq = 2(U — Uq). Moreover, the total energy 
E is related to the mean square radius of the atom cloud 
(r 2 ) according to E = 2E pot = Nmuj 2 (r 2 ). Since (if ) 
can be measured in a model-independent way [271 . [40| . 
the total energy E = E(T) as a function of the tempera- 
ture T provides an alternative method to determine the 
temperature of the interacting Fermi gas. We obtain a 
curve which is similar like Fig. 3J However, the major 
drawback of this method is that it requires knowledge of 
the ground state energy Eq, which must be determined 
with sufficient accuracy. As evident from Fig. [3J there 
are significant discrepancies in the ground state energy 
Eq for different experiments and theories. 
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VI. LOW-ENERGY COLLECTIVE MODES 

In the superfluid regime at very low temperatures the 
fermionic quasiparticles are frozen out and do not con- 
tribute to thermodynamic quantities because they have 
an energy gap which is related to the binding energy of 
the Cooper pairs. However, the spontaneous symmetry 
breaking implies a gapless Goldstonc mode which is the 
Bogoliubov- Anderson mode. This mode propagates with 
a constant velocity c like phonons. For this reason, the 
low temperature behavior of the entropy density and the 
internal energy density is ruled by the well known Stefan- 
Boltzmann formulas 



s(r) n(r) 



[u(r) -u (r)]n(r) 



8 g(r) T 3 
3 c(r) 

9 CT ( r ) ^4 

c(r) 



(6.1) 
(6.2) 



for phonons with one polarization degree of freedom 
where cr(r) = (7r 2 fc|)/(60fr 3 [c(r)] 2 ) is the Stefan-Boltz- 
mann factor. Since the sound velocity c(r) depends 
on the local particle density, er(r) and c(r) are space- 
dependent parameters. Eliminating the temperature we 
obtain a relation between the local entropy per particle 
s(r) and the local internal energy per particle u(r), which 
can be written in the form 



s(r) _ 2tt 
l^T ~~ 3 x 5 1 / 4 



vf{t) u(r) - M (r) 
c(r) e F (r) 



3/4 



(6.3) 



Here v F (r) = hkp{T)/m is the local Fermi velocity, and 
uq(t) is the ground-state energy per particle. The ratio 
c(r)/v F (r) = c/vp = {i/i) 1 ^ 2 = 0.345 is constant and 
related to the Bertsch parameter £ = 0.358. 

In our previous publication [24| we have calculated all 
thermodynamic quantities for the homogeneous system. 
Hence, also the function s = s(u) is available, so that 
we can check the asymptotic formula (|6.3p for the ho- 
mogeneous system. It turns out that the resulting ex- 
ponent at low energy is indeed equal to 3/4. However 
the amplitude, which is determined fully by the sound 
velocity c, gives c/v F = 0.7 for the unitary Fermi gas. 
This is about a factor of two larger than the expected 
value c/vp = (S./3) 1 / 2 from the Bertsch parameter. This 
discrepancy indicates, that the accuracy of our theory 
at very low temperatures is not sufficient to extract a 
reliable value of the sound velocity from the entropy. 

In the trapped case, unfortunately, asymptotic for- 
mulas like (|6.1|) - (|6.3|) do not hold for the total entropy 
S and the total energy E. To see this, wc integrate 
Eqs. 1|6.1[) and (|6. 2|) over the whole space. Using Eqs. 
(|2.5p and (|2.3p we obtain well defined results S and 
U — Uo = \(E — Eq) for the left hand sides, respectively. 
However, while the temperature T is constant, the space 
dependence on the right-hand sides arise from the factor 
<7(r)/c(r) ~ [c(r)]~ 3 ~ [v F (r)}- 3 ~ [n(r)]" 1 . This factor 
is minimum at the center of the trap but maximum at the 
surface of the atom cloud. Thus, the main contribution 




-2-10 1 
log 10 [(E-E )/NE F ] 

FIG. 5: (Color online) The local exponent b = d\nS/dln(E- 
Eq) of the function S(E) where E is two times the potential 
energy (green dashed), two times the internal energy (red dot- 
dashed), and one times the total energy (blue solid). Here Eo 
is the ground-state energy for T = and 5 = 0. 



of the integral arises from the surface of the atom cloud 
(see also Fig. [5]). Here the Fermi gas is a normal fluid, 
and the dimensionless temperature 6(r) = kBT/kp(r) is 
large, so that the Stcfan-Boltzmann formulas are not ap- 
plicable. Consequently, an asymptotic formula like (|6.3p 
cannot be derived for S and E — Eo of the whole trap. 

Empirically, it has been found in the experiments by 
Luo el al. [27[ that the total entropy at low energies varies 
with an effective power law 



S/Nk B ~ [(E - E )/NEp} 1 



(6.4) 



with an exponent b ks 0.59. In order to check whether 
such a behavior is consistent with a microscopic theory, 
we determine the exponent b = din S/dln(E—Eo) by log- 
arithmic differentiation of the blue-green-red solid curve 
S = S(E) in Fig. [3] The result is shown as blue solid 
curve in Fig. [5] and confirms that S(E) in a trap indeed 
exhibits a power law behavior in the regime near the 
ground state. We have adjusted the ground-state energy 
Eq by fine tuning in order to obtain a well defined ex- 
ponent in the limit E — > Eq. The resulting exponent 
b = 0.70 is surprisingly close to the value 0.75 of the lo- 
cal asymptotic formula (|6.3p but differs from the value 
b = 0.59 ± 0.03 inferred from the experimental fit. 

In Fig. [5] the superfluid transition is located at the 
position log 10 [(£ - E )/NE F ] = -0.65. The left part of 
the blue solid curve corresponds to the superfluid region. 
Here the exponent is nearly constant up to the superfluid 
transition. Even though the Stefan-Boltzmann formulas 
are not valid, the exponent b = 0.70 is remarkably close 
to the theoretical value 0.75. The right part of the blue 
solid curve corresponds to the normal fluid region. Here 
the logarithmic derivative of the entropy with respect to 
energy decreases monotonically with increasing energy E 
and thus no well defined exponent can be attributed to 
the high energy part of the curve. 
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The virial theorem implies the energy relations E = 
2-Bp 0t = 2U. In order to check the validity of the 
virial theorem we consider the entropy functions S = 
S(Ep t) and S = S(U) and calculate the related expo- 
nents b(E pot ) = d\nS/d\n(E pot - E pa t,o) and b(U) = 
d\nS/dln(U — Uq), which are shown in Fig. [5] as green 
dashed line and as red dot-dashed line, respectively. 
These lines should be compared with the blue solid line, 
which represents the exponent b(E) = d In S/d \n(E — 
Eq). In the normal- fluid region the virial theorem is well 
satisfied, because the right parts of the curves are lying 
nearly on top of each other where the small deviations 
are numerical errors. 

In the superfluid region the virial theorem is satisfied 
only approximately because of the modification of the 
theory in order to have a gapless Bogoliubov- Anderson 
mode (see Sec. II. J in Ref . 1241) . As a consequence, the left 
parts of the curves deviate from each other. We find three 
different values for the, exponents which are b(E pot ) = 
0.75, b(U) = 0.65, and b(E) = 0.70. These results are 
related to the three different ground-state energies (|3.4j) . 
p.5p . and (|3.6p . respectively. 

VII. CONCLUSIONS 

Based on our previous results for the BCS to BEC 
crossover problem in a homogeneous gas [24[, we have 
calculated density and entropy profiles in a trap within a 
local density approximation. In addition, we have deter- 
mined the total entropy S and energy E in the unitary 
regime and have compared our results with both exper- 
iment and recent theories. For temperatures above the 
superfluid transition temperature, very good agreement 
is obtained. However, the value of the critical temper- 
ature and the behavior at very low temperatures differ 
appreciably from the experimental estimates and their 
theoretical analysis in earlier work. 

First of all, our value for the Bertsch parameter £ and 
thus the ground-state energy Eq is about 10% smaller 
than the results obtained from variational Monte Carlo 
calculations or from Gaussian fluctuation theories around 
the BCS ansatz for the ground state. While our value 
£ = 0.36 agrees well with the most precise results so 
far obtained from an e = 4 — d expansion [35| . it dif- 
fers from those obtained from variational Monte Carlo 
calculations [13, [H[, or from those including Gaussian 
fluctuations around the BCS mean field [28l l39j. where 
£ = 0.42(1) or £ = 0.40, respectively. Given the uncer- 
tainty in determining £ experimentally (which requires an 
extrapolation to zero temperature) it is clearly important 
for theory to make precise predictions for £ which do not 



rely on approximations that are apparently limiting all 
present results. In view of the fundamental importance 
of this parameter in the context of strongly interacting 
Fermi gases, progress here would be highly desirable. 

As a second point, our values for both the critical tem- 
perature and value of entropy at T c are appreciably lower 
than those obtained in the theories of Bulgac et al. [3(| 
and of Hu et al. [29[ and also those inferred from the 
original analysis of the experimental data [27|. Now, as 
is evident from Fig. [31 a measurement of the function 
S(E) does not provide a sensitive measure of the critical 
temperature. Quantitative results for the critical tem- 
perature of the unitary gas have been obtained by Shin 
et al. (53| . They rely on using gases with a finite imbal- 
ance nj 7^ which have always a fully polarized outer 
shell in a trap. Since a single species ultracold Fermi 
gas is nonintcracting, its temperature can be reliably de- 
termined from cloud profiles. An extrapolation back to 
zero imbalance gives a critical temperature at unitarity 
which is close to the value predicted both in our the- 
ory and in Monte-Carlo calculations by Burovski et al. 
[HI, EH and Bulgac et al. [48| . As pointed out in section 
IV, there is now evidence that our result fc^T^ ~ 0.21 Ep 
for the critical temperature of the unitary gas in a trap 
is rather precise. Together with the entropy-temperature 
curve shown in Fig.[2J this would allow doing precise ther- 
mometry for balanced gases, that has not been possible 
so far. 

Finally, we have shown that for temperatures of order 
fcgT w 0.06 £f(0) ~ 0.10 Ep which have been reached 
experimentally [53j], an efficient way of further cooling 
the gas is possible by removing the high entropy outer 
part of the atomic cloud and readjusting the confining 
potential. This method is similar in spirit than the stan- 
dard evaporative cooling idea, but potentially much more 
efficient. It might open the avenue to reach regimes in 
which the entropy per particle is much less than ks, a 
necessary condition for realizing many of the nontrivially 
ordered states that are in principle accessible with ultra- 
cold fermions [22| . 
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